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We expand the solutions of linearly coupled Mathieu equations in terms of infi- 
nite continued matrix inversions, and use it to find the modes which diagonalize the 
dynamical problem. This allows obtaining explicitly the ('Floquet-Lyapunov') trans- 
formation to coordinates in which the motion is that of decoupled linear oscillators. 
We use this transformation to solve the Heisenberg equations of the corresponding 
quantum-mechanical problem, and find the quantum wavefunctions for stable oscilla- 
tions, expressed in configuration-space. The obtained transformation and quantum 
solutions can be applied to more general linear systems with periodic coefficients 
(coupled Hill equations, periodically driven parametric oscillators), and to nonlinear 
systems as a starting point for convenient perturbative treatment of the nonlinearity. 
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I. INTRODUCTION 

The Mathieu equation for a single degree of freedom is very well known [l[ . In this paper 
we discuss the coupled system of Mathieu equations 

ii+[A-2Qcos2t]u = 0, (1) 

where u is an /-component vector and A and Q are constant symmetric / x / matrices. We 
also extend the treatment of eq. ([T]) to include an inhomogeneous right-hand-side, and more 
general vr-periodic coupled parametric oscillators. 

This multidimensional matrix equation has been researched thoroughly (see e.g. ^ and 
[sl where also some applications are exemplified, and 0^)- Many general treatments of this 
system are perturbative and concerned with stability analysis - i.e. with finding the regions 
in some parameter-space for which the solutions are bounded. In this contribution we are 
interested primarily with describing the classical and quantum solutions of the system in 
terms of decoupled modes of oscillation. The solutions presented here may find application 
in the description of various discrete systems of coupled parametric oscillators, e.g., among 
others, trapped ion crystals 0, S], coupled arrays of nanoelectro mechanical oscillators [q, 10 
and binary Bose- Einstein condensates 111 . 

The quantum problem of time- dependent linear and quadratic Hamiltonians has also 



been considered in many publications. One of the first treatments is by Husimi [12], who 
considered the problem of the one- dimensional quantum parametric oscillator with a driving 
force. He constructed Gaussian wavefunctions assuming that the classical solutions are 
known, and obtained the propagator of the system. He also derived transition amplitudes 
between states of Hamiltonians which are time-independent at some initial and final times. 

The one- dimensional problem has been especially important to the description of the 
motion of single ions trapped in radiofrequency traps. It has been treated extensively using 



different methods and summarized in a few texts, e.g. 13|, and see references within. 



In [15[ Lewis and Reiesenfeld have considered a general time-dependent Hamiltonian, 
and have shown that the eigenvalues of an invariant operator (whose total time-derivative 
vanishes), are time- independent, and that the eigenstates can be chosen with a specific time- 
dependent phase, so as to solve the Schrodinger equation. This theory is the basis for most 
treatments of multidimensional time-dependent Hamiltonians. 



In 



161] coherent states and eigenfunctions have been obtained for a diagonal system of 



harmonic oscillators with a time-dependent frequency. Holz [17| has considered multidi- 
mensional time-dependent oscillator Hamiltonians which remain positive-definite. He has 
constructed the coherent states, assuming that the Lewis-Riesenfeld invariants are known 
and that the Hamiltonian at the initial time is time-independent. Transition amplitudes 
have also been calculated in these works. In [isj expressions for the wavefunctions of gen- 
eral time-dependent upto-quadratic Hamiltonians are given formally using solutions of the 
classical equations in phase-space. Leach [l9| has considered the use of a time-dependent 
transformation to obtain a time-independent Hamiltonian, for the case of positive-definite 
Hamiltonians, and in terms of formal classical solutions. In this contribution we give the 
coherent and number states of the Schrodinger equation explicitly in terms of the decoupled 
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modes of eq. ([Tj). 

This paper is organized as follows. In Sec. [Ill we obtain an analytic expansion for the 
solutions of eq. ([1]), which is not based on a small parameter, but rather uses infinite contin- 
ued matrix inversions. We obtain explicitly a time-dependent transformation to coordinates 
in which the motion is that of decoupled linear oscillators. In Sec. IIIII we use this trans- 
formation to find the wavefunctions of the corresponding quantum system. We conclude 
in Sec. HV] with a summary of our results and comments of possible further research and 
applications to various physical systems. The classical linearized modes of an ion crystal in 
a Paul trap are treated in |20| using the methods described here. 

II. SOLUTION OF THE LINEAR EQUATIONS 
A. The Floquet Problem 

Eq. (Il]) is a homogenous linear differential equation with periodic coefficients and therefore 
amenable to treatment using Floquet theory. In this subsection we recall a few facts from 
this theory and introduce the Floquet-Lyapunov transformation which will allow us to obtain 
explicitly the classical and quantum solutions. 

Eq. ([T]), with a vr-periodic parametric drive, can be obtained by a suitable nondimen- 
sionalization (scaling) of the coordinates and time of the equations of motion (e.o.m) of a 
physical system parametrically driven at frequency Q. For the Newtonian problem with / 
degrees of freedom, the corresponding Floquet problem is stated in terms of coordinates in 
2/-dimensional phase space by defining 

where 1/ is the /-dimensional identity matrix. The e.o.m is written in standard form as 

= n(t)0. (3) 

In the following, an /-dimensional vector u will be denoted by a lower case Latin letter 
(usually with an arrow), and Latin subscripts will be used for its components (um)- f- 
dimensional matrices will be denoted by capital Latin letters (Q). A 2/-dimensional vector 
will be denoted by a lower case Greek letter (with no arrow), and Greek subscripts will 
be used for components of 2 /-dimensional vectors (0^). Capital Greek letters (unitalicized) 
will denote 2/-dimensional matrices (11, or B). 

A fundamental matrix solution to eq. ([3]) has 2/ linearly independent column solutions 
and obeys the matrix equation 

$ (t) = n (t) $ it) . (4) 

A principal fundamental matrix solution $ (t) is a fundamental matrix solution which equals 
the identity matrix at some point in time. We will always take a principal fundamental 
matrix solution to obey this at t = 0, i.e. 



$ (0) = l2f. 



(5) 
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and then $, which is obviously unique, is also known as the matrizant of eq. (|3])). 

Now let T = 27r/fi = vr be the period of H (t), i.e. U{t + T) = U{t). Then $ (t + T) 
is also a fundamental matrix solution. Therefore its columns must be linear combinations 
of the columns of $ (t), i.e. $ (t + T) = $ (t) S for some non-singular constant matrix S. 
In particular, given the initial conditions in eq. ([5]), we find that H = $ (T), known as the 
monodromy matrix. H can be brought to Jordan canonical form, which we assume to be 
diagonal (which holds for the case of stable oscillations). The diagonalization is given by 

p-iSP = A = diag{Ai,...,A2/} (6) 

where are the complex Floquet eigenvalues (also known as multipliers). 
Applying the coordinate transformation by writing 

$ (t) = T (t) P^^ (7) 

we get that T (t + T) = T {t) P^^SP = T it) A. Therefore the column-vector of T [t] 
obeys {t + T) = XyTy (t), and thus T^, (t + riT) = A"T,^ (t). Consequently the solutions 
are decaying for \X^\ < 1 and unstable for \Xu\ > 1. When eqs. (E])-© are Hamiltonian 
(as in our case), the sets {X^} and {(A*)^ } must coincide. Since the equations have real 
coefficients, nonreal X^, must come in conjugate pairs. Therefore nonreal X^ not on the unit 
circle come in quadruplets as X^, (A*)~ , A*, A~^ and real Aj, not on the unit circle come in 
pairs as A^,, A~^. Letting Aj, = e*''"^ defines the characteristic exponents = -^InX^. We 
multiply Tj, (t + T) = X^T,^ (t) by e"*^"*-*"*""^^ and obtain the normal form for the column- 
vector solutions 

(t) = (t) e^^"\ (8) 

where are T-periodic vectors, F^ (t + T) = T^, (t). 
Using eqs. (I7j) and (E]) and defining 

B = diag{z/3i,...,^/32/}, (9) 

we may now write 

$ (t) = F [t) e^*p-^ = F (t) e^*F-i (0) . (10) 

We note that in the general case, P^^ = F^^ (0) will not be unitary, although this may happen 
in certain highly symmetric cases, e.g., in the trivial case in which A and Q commute and 
there exists a constant orthogonal transformation which diagonalizes eq. ([I]) into a system 
of decoupled Mathieu equations. 

Differentiating eq. ( ITOj) and substituting into the e.o.m, eq. (jl]), we have 

(te^t + F5e^*) P-^ = nFe^*p-^ (11) 

or, 

f + FB = HF. (12) 

If we now substitute into eq. the time-dependent coordinate change 



<P{t)=T{t)x{t), 



(13) 
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we get 

tx + Tx = UTx, (14) 
which, after using eq. ( !T2|) and multiplying on the left by T~^, reduces to 

X = Bx. (15) 

Eq. ( |T5l) is a diagonal equation with constant coefficients, the solutions of which are the 
Floquet modes 

Xu (t) = X. (0) e^^^*. (16) 

The time-dependent transformation T (t) of eq. ( lT3l) is known as the Floquet-Lyapunov 
transformation. 



B. Solution Using an Expansion in Infinite Continued Matrix Inversions 

We turn to an analytical expansion of the solutions of the homogenous e.o.m ([1]). The 
following expansion allows to obtain the frequencies and the coefficients of the solution 
vectors in a generalization of an infinite continued fractions expansion, to arbitrary precision. 
Infinite recurrence relations have been used for solving various types of differential equations 



(see e.g. chapter 9 of |2l|) and differential- delay equations |22L and applied recently to the 
study of the stability of a trapped Bose- Einstein condensate [23|. The method described 
below gives the solution in a form which is immediately suitable for obtaining the Floquet- 
Lyapunov transformation. 

We seek for solutions of eq. ([T]), in the form of a sum of two linearly independent complex 
solutions, 



u 



where b and c are complex constants determined by the initial conditions. In general, 
the characteristic exponents P may be complex. Except when there are /3's which are 
real and integral, a system composed of solutions of the form of eq. (fT7|) will constitute a 
fundamental system. For an integral (3, eq. (fT7|) gives a single vr- or 2'7r-periodic solution and 
the other linearly independent solution will in general (excluding the trivial case Q = 0) 
be unbound [l|, We do not treat this marginal case as we will soon restrict ourselves to 
stable oscillations only. Following Sec. Ill A| stable modes will be described by f3 taking a 
real nonintegral value, which can obviously be chosen in the range < P < 2, /3 ^ 1. For 
the stable modes, when /3 is real, we find that b = c* and are all real. 

We assign eq. ( IT71) in eq. ([T]), discard the negative exponent terms (which give identical 
relations), and find that the solutions must obey for all t 

- J2 C2n {2n + e*(2n+/3)t + _ g (^^2* ^ g-,,2i)] J2 C72ne^('"+'')* = 0, (18) 

where in the above expression and for the rest of this section, the summation is over n G Z. 
Thus we get the recursion relation 

- C2n (2n + Pf + AC2n - Q (^2„-2 + (?2n+2) = 0. (19) 
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By defining 

R2n = A- {2n + , (20) 

we can write the first infinite recursion relation, which progresses towards positive values of 

n, 

QC2n-2 = R2nC2n ~ QC2n+2- (21) 

To get a second relation which progresses towards negative values of n, we reorganize eq. (j2Tll . 
obtaining 

QC2n+2 = R2nC2n ~ QC2n-2- (22) 

From these relations we can now obtain an expansion in infinite continued matrix inver- 
sions. Starting with n = 1 in eq. fl2T]) we get by repeated substitutions 

C2 = T2,^QCo = (jR2-Q[R^-Q [Re - ...r' Q] ' q] QCo. (23) 

Substituting decreasing values of n starting with n = in eq. ( 122|) . we get the independent 
relation 

QC2 = RoCo - QC^2 = f,,A =[R^-Q [R-2 - Q [R-A - Q] ' Q) Co- (24) 
Multiplying eq. ( |23|) by Q and defining 

Y2,^ = - QT2,pQ, (25) 

we find that all characteristic exponents P are zeros of the determinant of 1^2,/? (which is 
a function of /3). If there are degenerate /3's they will appear as degenerate zeros of this 
determinant. The vector Co for each /3 is an eigenvector of ¥2^13 with eigenvalue 0. Since A 
and Q are symmetric, l2,/3 is symmetric as well, and so its kernel will be of dimension equal 
to the algebraic multiplicity of the /3 root. The vector C2 can be obtained by an application 
of T2^i3Q to Co, for 72 = —1 we use C_2 = [T~2,/3]~ QCq, and so on for the other vectors. We 
note that the different vectors C2n,i3 are not orthogonal in general, and the vectors at every 
order in n mix different coordinates. 

Because of the presence of the diagonal term (2n + /3)^ in i?2n, we would have ||-R2nili 
(2n + /3) +0 (1), and therefore the general term of the expansion vanishes. Either A or 
Q may be singular and the expansion can still be applied in general. Even if both are 
singular, the expansion is valid if there are no integral values of (3, a case which we do 
not tackle as noted above. Excluding perhaps isolated values of f3 (and atypically in the 
parameter-space), all matrices which are inverted in the above expressions will be invertible, 
and while employing the algorithm in practice, the invertibility of the matrices is of course 
easily verified at each step. In App. [A] we extend the infinite matrix inversions to obtain 
the periodic solution of eq. ([1]) with an inhomogenous r.h.s, and also comment on some 
computational aspects of this method. In App. [B]we show briefly how the method may be 
extended to a system of coupled Hill equations. 
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C. The Floquet-Lyapunov Transformation For Stable Modes 

We can now find explicitly the time-dependent Floquet-Lyapunov transformation of 
eq. f lT3|) which transforms the Floquet problem to a time-independent equation, which in 
our construction is also diagonal. We further assume that all Floquet modes are stable, i.e. 
that the 2/ linearly-independent solutions of eq. ([3]) are oscillatory and thus characteristic 
exponent come in complex conjugate pairs. This simplifies many expressions and avoids 
complications in the quantization, since the eigenf unctions of the negative harmonic poten- 



tial (the parabola potential 2^) are not square integrable over the real line. We therefore 



take B of eq. ([9]) in the block form 

^={'o -w)' (i?)^,^ = diag{A,...,/3^} (26) 

where Pj are positive. We define the /-dimensional matrix U whose columns are constructed 
from the series of /-dimensional vectors C2n,i3j obtained from the recursion relations for the 
solutions u of eq. ( |T71) . i.e. 

(^)/x/=(E(?2„,,.e^^"* ...) (27) 

We similarly define the /-dimensional matrix V composed from column- vectors as 

(^)/x/ ={^E{2n + P,) C72.,;,,e^2«* ... ) , (28) 

and thus we may represent a 2/-dimensional fundamental matrix solution in the form 
\E' (t) e^* where \E' (t) is written in block form 

where U* denotes the complex conjugate (and not hermitian conjugate) of the matrix U. 

By multiplying (t) e^* on the right with \E'~^ (0) we get the matrizant $ (t) since the 
initial condition of eq. ([5]) is obeyed. Then by comparing with eq. f lTOj) we find 

^ (t) e^*^-i (0) = r (t) e^*r-i (0) , (30) 

so that we may choose 

r(t) = *(t)=ff! ^l) (31) 



(the choice is in fact unique only up to a constant matrix which commutes with B, and we 
will use this fact in the following). 



III. QUANTIZATION 
A. Hamiltonian Formalism 



In this subsection we consider the results of the previous section within the Hamiltonian 
formalism. We find the conditions for the Floquet-Lyapunov transformation to be canonical. 
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which also allows us to obtain its inverse explicitly in terms of matrix transpositions and 
complex conjugation. This requires finding the generating function of classical mechanics, 
and gives the transformed Hamiltonian. These results are novel to the best of our knowledge. 
Let us rewrite the e.o.m, eq. ([3]), in the form 

J0 = jn0 = H0, (32) 

where 

f A-2Qcos2t 

1/ J' 1, 

J is a skew-symmetric matrix (J^^ = J* = —J) and H is symmetric. Denoting hj p = u the 
momenta canonically conjugate to the coordinates u, eq. (l32l) is seen to be the canonical 
form of Hamilton's equation, 

with the corresponding Hamiltonian "H written as the quadratic form 

n = ^0*H0, (34) 

where denotes the transposed vector. 

Similarly, we rewrite the transformed e.o.m, eq. f|T5|) . in the form 



Kx = KBx, (35) 
where K is the antihermitian matrix (K^^ = = — K) 



K 



-ilf 
^lf 



and KB is a hermitian matrix (in fact, positive definite), with the explicit form KB = 
diag{/3i, .. .,/?/}. 

Let us here introduce explicitly the connonically conjugate variables of the Hamiltonian 
formalism (we here break the notation a little). 



X 



where ^ is the new /-dimensional vector of coordinates and we see that —i( are the new 
conjugate momenta, such that eq. fl35l) is derivable from the Hamiltonian 

n'^^x'i^x + y{t), (36) 

with 
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B is given by eq. ( 126|) and y (t) is a function of time alone, that does not enter the equations 
of motion. In the following we will prove that indeed H' is the transformed Hamiltonian 
provided that the Floquet-Lyapunov transformation is canonical. Our choice of the Floquet- 
Lyapunov transformation will be such that classically, y (t) = 0. Using the realness of u 
and p and the explicit expressions for T (t) and (t) it is easy to verify that 

e = C- (38) 

In App. [O we obtain an expression for the inverse of the Floquet-Lyapunov transforma- 
tion, and show that the matrices U and V can be rescaled by multiplication with a (diagonal) 
matrix, such that 

U {-2iV' (0) U (0)) , (39) 
and V accordingly, thereby imposing the following normalization condition, 

V (0) U (0) = \i. (40) 
By eqs. ([02])- ([05]) subject to eq. (gO]), we get 

and by writing F'^F = I2/ in block form, we find the identities 

jjty* _ ytu* = jjty ^ ytjj (^42) 

which we will use below. 

We now turn to finding the classical generating function of the canonical transformation 
relating T-L to %' . If x (t) = T"^ (t) (t) is to be a canonical transformation, we search 
for the generating function {u, t), expressed in terms of the old coordinates and new 
momenta, which obeys the set of equations 

and the transformed Hamiltonian is given by 

H' = H + dJjdt. (44) 
With the help of eq. fHT]) we can invert for p in terms of u and C to get 

p = -lU-'C + U'^Vu, i = iV^u - iU^ {U-'V'u - lU-X) , 

where we define for bravity f/^* = [f/^^]*, and in the following we will use and 0~* with 
a similar definition. A solution to eqs. ( H3|) exists and reads 

jr = lu*f/-V*M - lu'U'X + ]^iC,'U^U-\ (45) 

provided that 

- U^U-*V* = -iU-^ (46) 

and that U~'^V^ and WU~^ are symmetric matrices. These conditions follow after some 
manipulations from the identities in eq. fB2|) . Thus we see that the normalization of eq. f HOj) 
guarantees that the Floquet-Lyapunov transformation is canonical. 
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B. Quantization 

In this subsection we apply the Floquet-Lyapunov transformation to the operators in the 
Heisenberg picture of the quantum problem corresponding to eq. ([1]), allowing to diagonalize 
it in terms of ladder operators. We then find explicit expressions for the wavefunctions of 
the coherent and number states, utilizing the periodicity of the Floquet-Lyapunov trans- 
formation and its inverse, and the decoupled oscillatory modes with their characteristic 
frequencies. These results are novel to the best of our knowledge. 

We first canonically quantize the system by promoting the canonically conjugate variables 
u and p to operators obeying the quantum commutation relations 

[Pj.Uk] = -iSjk, [Pj,Pk] = [uj,Uk] = 0, (47) 

where we will denote operators with a hat, and set h = 1. The Heisenberg equations of 
motion for these operators, 0, are identical to eq. (j32l) . Repeating the derivation of Sec. Ill Al 
we see that the noncommutativity of the operators has no affect on the transformation and 
thus we find in the Heisenberg picture the e.o.m 

= Hx, (48) 

which, using eq. ( 126|) . is the diagonal set 

4- = ^Pji 4 = -^/3.C„ (49) 

the solution of which is simply 

4- (t) = 4- (0) e^/'^*, C,(t)=C,(0)e-''^*. (50) 

By assigning eq. (HTj) in eq. (H7I) . we get that the canonical commutation relations of 
these operators obey 

Cj (t) , 4 (t)] = (51) 

with all other commutators zero, and this result is subject to the normalization condition 
eq. (HOl) which ensures that the Floquet-Lyapunov transformation is canonical. 

The commutation relations of eq. (15T|) are easily recognizable as those of the creation and 
annihilation operators of a harmonic oscillator, since the hermiticity of u and p immediately 
implies = (j (which also follows directly from eq. ( [38|) ). We may therefore define the 
time- independent eigenstates of ( (t), the coherent states, in the Heisenberg picture, by 

0(^)10 = 0(^)10 = 0(0)6-^/^^*10, (52) 

and the normalization and completeness relations are 

(CIO)=e^*<', /rf/i/IO(0 = i, (53) 



where dfif = n ^ e ^*'^d^ and C, = Yi dxj Y[ dyj if C is written in terms of the real variables 

~ "''i + "^Vj ■ 
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We next show that the Schrodinger wavefunction of the coherent state vector of the 
system, is given in the coordinate representation by 

(m, t) = {u\C) = U exp {u, C, } , (54) 

where J-" {u, ( (t) , t) has the same functional form of the classical generating function 
J-'{u,(,t), only with the explicit time-dependence of ( = C{t), which we will omit be- 
low, except where necessary. Multiplying eq. (152|) by {u\ we obtain the differential equation 



-2 J2 + Yl Ujkdu, ) {u I c) = {u I c) • 

k k / 



(55) 



Since matrix elements are invariant under any unitary transformation, eq. (155|) holds in the 
Schrodinger picture as well. The solution of this equation is 

{u\0=Af exp I hu'U-'V'u + u'U-X + / (C) | • (56) 

To determine the (time-dependent) normalization constant A/" and the function / {(), we 
first impose the normalization of eq. ( 153|) . 



^"^' = (ClC')= / {C\u){u\C)d^u. (57) 



To evaluate eq. ( 1571) we use the following result. Let T be a symmetric n x n complex 
matrix with positive definite real part, and b a complex vector. Then 

J = j exp|-ix*Ta; + 6*x|d"x = (27r)"/^(detT)-^/^exp|i6*T-i6|, (58) 

and the value of (det T)~^^^ is defined by analytic continuation, writing T = 9^eT + ieJmT 
starting with the (positive definite) real part of T and increasing e continously to 1. 
From eq. fj42|) we get 

f/-V* = U-^V^ + tU-^U-\ (59) 



so 



e^*-^' = £^ exp |-ln*t/-tf/-iu + u^U^'C + u^U'^C + f (C) + /* (C)} (60) 
and we get (since {U'^U) ^ is obviously positive definite). 



e<*<' = i2ny/' det {UU^Y^' 



X 



X 



exp |i {Cu^u-'C + c'u-*uc + C'U'*uu^u-\' + ec) + / (CO + r (C)} • (61) 
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By using the fact that WU * is symmetric we find 

(U-*Uy = U^U-\ U-*UU^U-' = If, 

so we may deduce 

/(C) = -^C'U^U-X, \^f\ = (27r)--^/^det (f/f/t)'^/^ 



(62) 



(63) 



which proves eq. fl54|) . without fixing the phase of Af. 

Finally, to determine the wavefunction completely, we require that {u, t) be a solution 
of the Schrodinger equation, which is expressed in coordindate u space using the Hamiltonian 
of eq. (El, 

z^tV^c (^' t) = ^V'c (w. t) ■ (64) 



Substituting eq. (iMll in eq. ( iMll and inserting a resolution of the identity from eq. (]53ll we 
have 



j^/Ar + idtT) {u\o 



d^'f{u\C) {C\n\c). 



In App. Owe show that the integral on the r.h.s of eq. ( I65l) is equal to 



so that we get 



-dt^+^ti{B-iW} {u\C) 



Af/Af=-^tT {iB + W}. 



(65) 



(66) 



(67) 



Writing Af = \Af \ exp {i arg A/"}, the above equation is equivalent to the two equations 



Af 



/\^f \ = --^e{tTW} 



and 



dtHTgAf = -- {tiB + 3m {trW}) 



(68) 



(69) 



Eq. fl68|) is in fact an identity which results from eq. f l63|) as shown in App. [Cl where it is 
also shown that the solution of eq. (|69|) is 



arg 



Ar=-i^/3,t-iargdetf/. 



(70) 



Let us write here again the complete expression for the coherent state vector if(^, 



= (27r)~-^/^det (f/f/t) 



-1/4 



X 



X exp 



(71) 
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where ( is the vector with components Q (0) e~''^^*. 

In the case of a single harmonic oscillator of frequency (3, eq. (!7T|) must reduce to the 
familiar wavefunction in configuration space of a coherent state, with complex label (q = 
( (0). This can be seen by noting that for this case, the transformation matrices U and V 

with the normalization of eq. ( I40l) . become the scalars U = and V = 'i\f^, and then 



eq. (I7T1) becomes, with m = 1, ^ = 1, the expected expression (e.g. eq. (21.1.132) of |25| ) 



A = (2vr) 



-1/4 



2/3 



-1/4 



exp 



^f3t - ^f3u' + y^uCoe-'^' - ^C^^' 



We will now find a complete orthonormal basis of solutions of the Schrodinger equation. 
For that purpose we can use a generating function for multidimensional Hermite polynomials 
H9 , defined by 



Gc = exp <! x'C( - -CCC 



^ fill nA 



(72) 



where C is a symmetric matrix, and the summation is over all /-tuples of nonnegative 
integers, n. An explicit definition of H9 is given by 26 



H? (x) 



-x*Cx/2 



Then, setting in eq. fl72|) x = U'*u and 

c = u^u-\ 

we can write ip^^ as 

.Ci (0)"^ 



(73) 



(74) 



1 



V'c = A/'exp <^ 2 V u} > e 



ni! 



C/ (0)^ 



-if? {U~*u) . (75) 



Eq. fl75|) can be interpreted as an expansion of ipt^ in terms of a complete orthonormal set of 
solutions of the Schrodinger equation, ip^, and in that case the coefficients of the expansion 
must be time- independent. We may therefore write 

V-n = ^fne-'^^ "^''^^ exp I hu'U-'V'u^ {U-*u) , (76) 

where Mh = CfiM and is time- independent, and we impose the normalization 



In App. Owe show that 
so that we have 



[ni\ ■ ■ - nf. 



-1/2 



(77) 



(78) 



14 



(2n) 




X exp 



det [UU^y 



-1/4 



X 



|-2 ^ + - arg det [/ + ^m^^y-^y^wj {U' 



u 



■). (79) 



Let us note how to get from eq. (ff9l) the famihar exp ression for the wavefunctions of the 



one- dimensional case, e.g. in the form of eq. (36) of [13|]. The latter is given in terms of the 
periodic function $ {t) and the constant v defined there, in eqs. (22)-(23), with {3 = (3^. Let 
us keep for now the nondimensional units (with the drive frequency being equal to 2), so 
z/ = ^ (2n + (3) C2n- By the normalization of eq. fHOll . we see that U (t) = (t). The 

usual one- dimensional Hermite polynomials Hn are obtained in eq. (!73|) by setting C = 2. 
Since we have C = $*/$, this requires the variable change 



which gives that 



In addition we have 




d_ 

dx 




2-l/2g-iarg$ 



-in are <I> 



(27r)-'/^ (det UU^) 



d_ 
dz'' 



-Hr, 



-1/4 g-liargdetC/ 



$1 



tU 



1/4 



$1/2 



and ^iU~^V^ = ^§ ~ Eq. ([T]) is expressed in terms of rescaled time, and we now return 
to the time variable before the Mathieu scaling t — ?• Qt/2 (with Q being the physical drive 
frequency), and therefore put u — )■ z/fi/2, and put back h (we still have m = 1), to get the 
wavef unction 27 



g-marg* ^ 1/4 f / 1 \ ^ 1^ 1 A $ n ^\ 2] rr f \ 

= r + 2 J 2 ' + 2n 1^^$ - ^^-2 j " I V 

(80) 

In a multidimensional problem, two distinct situations may arise. If U is diagonal, which 
means that C is diagonal as well, the generating function of eq. ( 1721) obviously factorizes 
into a product of exponents, and the wavefunction will be a product of one-dimensional 
wavefunctions, each depending exclusively on one variable, as obtained above. As mentioned 
in Sec. Ill At U can be made be diagonal if there exists a constant matrix which diagonalizes 
the equations of motion. Then U will be diagonal in some normal modes which are time- 
independent linear combinations of the original coordinates. If such a diagonalization does 
not exist, the wavefunctions will depend on (time-dependent) complex linear combinations 
of the coordinates, through the multidimensional Hermite polynomials. 
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C. The Inhomogenous Equations 



In this subsection we describe briefly how to obtain the wavefunctions of the quantum 
system which corresponds to eq. ([1]) with a driven r.h.s in the form 



[A-2Q cos2t]u = G + 2Fcos2t, 



(81) 



where G and F are /-component constant vectors. We rewrite eq. (IHTi) in Floquet form 
(eq. (EJ) using 



G + 2F cos2t 



A 



as 



(j)-U{t)(f) = X{t) , (82) 

which we transform using the Floquet-Lyapunov transformation eq. ( ITSll and its inverse, to 
get the e.o.m in the form 

X-Bx = T (t)-' A. (83) 

Each of the equations, eq. ( 183|) and eq. (IHT]) . has a unique vr-periodic solution (see App. [X|) . 
and these solutions are related by the Floquet-Lyapunov transformation. Eq. ( 183|) can 
be solved immediately term-by-term, however in App. |A] we use infinite continued matrix 
inversions to find directly the periodic solution of eq. ( IHT]) . in the form 



2ne 



i2nt 



To get the wavefunctions, we use the method of 
ordindate u space is now 



12 



^u^{A-2Qcos2t)u- [G + 2Fcos2t) -u 



(84) 

The Schrodinger equation in co- 
^. (85) 



Performing the (unitary) coordinate change x = u — we have dt — ?■ —Utt + dt, and by 
introducing the additional unitary transofmration. 



ijj = exp < ■ X + iUt, (t) > cp 



(86) 



the Schrodinger equation for ip becomes that of the nondriven problem, whose solutions are 
given in eq. fITT]) and eq. ( !79|) . The phase a-,, (t) is in fact the classical action of the vr-periodic 
solution, 



1 / • N 2 



t r 



-ui{A-2Qcos2t)u^ - (G + 2Fcos2t) ■ 



(87) 



which may be written more compactly using eq. ( l8T]l . 



«7r {t) 



*1 

2 



(ii^ + G + 2Fcos2t 



(88) 
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A closed algebraic expression for (t) can in fact be obtained by expanding the integrand 
of eq. (l88|l using eq. ( |8^ into a sum of exponentials, for which the integration is immediate. 

Thus, with (x, t) a solution of the nondriven Schrodinger equation, the solutions of 
eq. dHS]) will be 

ijj = exp < ■ {u — Ut,) + zajr [t] > (p[u — u^, t) . (89) 



IV. CONCLUDING COMMENTS 



Eq. ([T]) describes a coupled system of Mathieu equations. This equation can be considered 
as consisting of the first two terms in the expansion in a Fourier series of a more general 
system of coupled Hill equations [l|. 



u 



A 



n=l 



COS 2nt 



u 



(90) 



where the sum may be infinite in principle, and here we take the equation to be time-reversal 
invariant. The technique of expansion of Sec. IIIBI in matrix inversions can be applied 

The Floquet-Lyapunov 
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to solve eq. (]90ll with only the algebraic overhead growing 
transformation and the entire quantum treatment remain identical. In App. [B] we solve a 
Hill system with two harmonics. 

Beside the above generalization, we have considered in App. |A] an inhomogenous system 
(of a specific form) and its quantum counterpart is considered in Sec. IIII Ci Other types of 
driving can be handled similarly, and simple known transformations j3] can be used to handle 
more general linear system similar to eq. (190|) . such as systems with first-order derivatives 
(e.g. linear damping, gyroscopic forces or magnetic fields). 

Finally, we note that a linearization of a general nonlinear multidimensional system with 
periodic coefficients will lead to similar equations. After obtaining the Floquet-Lyapunov 
transformation, the system of decoupled time-independent oscillators can be canonically 
transformed to action-angle coordinates, for example. The obtained modes can be used as 
the zeroth-order approximation in a perturbative treatment of the nonlinearity 28|, l29 . 

An application of the methods described here to the analysis of the classical linearized 
modes of an ion crystal in a Paul trap can be found in 20| , where the accuracy of the solution 
is demonstrated by comparing to exact numerical simulations of the nonlinear problem. 
Various other linear and nonlinear parametrically driven physical systems can be accurately 
described and analyzed, either in the classical regime or as they are cooled close to the 
mantum ground-state of motion. Coupled arrays of nanoelectromechanical oscillators ^, 



lOl . ISOj are one exaraple. Parametric driving has also been recently applied to Bose-Einstein 



condensates [Ul, l3l|, l32| , and in particular the perturbations (i.e. phonons) of two-component 



condensates obey coupled Mathieu equations [33 
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Appendix A: Further Comments on The Infinite Continued Matrix Inversions 

Adding to eq. an inhomogeneous r.h.s we put it in the form 

{i+[A-2Qcos2t]u = G + 2Fcos2t, (Al) 

where G and F are /-component constant vectors. This equation will have a unique vr- 
periodic solution under some conditions [2] , a sufficient condition being that the homogenous 
equation does not have any vr-periodic solution (except the trivial one). This condition is of 
course fulfilled when the homogenous system has purely oscillatory modes. We first find a 



particular solution of eq. (lAip . using the method of Sec. Ill B[ We assign = J2nez -S2ne*^"* 



in the e.o.m and get 

{A - An^) B2n - Q {B2n^2 + ^a^+s) = G5nfl + F {5n,l + 5n,-l) . (A2) 

We write, defining = A — (2n)^ and using i?2n = -B-2n, 

ABo - 2QB2 = G (A3) 

R2B2 - g (So + = F (A4) 

R2nB2n - Q (^a^-s + ^Sn+a) =0, n>2. (A5) 

Eq. flASP immediately gives a recursion relation in the form of eq. fl2T]) (only i?2n here is 
defined differently), which allows us to obtain the expression in infinite inversions 

S4 = T2QB2, (A6) 

where 

T2 =[r,-Q [Re -Q[Rs- ...r' Q] ' Q] ' ■ (A7) 
Substituting eq. ( 1A6I) in eqs. ( lA3[) -( lA4l) we get the linear system 



A -2Q \ Bo] ^ G 
-Q R2 - QT2Q J \B2j \F 

which is readily solved, and the rest of the coefficients follow immediately. 



(A8) 
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We conclude by commenting on the computational aspects of the above method. The 
complexity of matrix multiplications and inversion, and computing a determinant are all 
equal, O (/^) or a bit better with improved algorithms. Using the method of matrix inver- 
sions, finding all zeros of the determinant of eq. (125|) can be done with complexity indepen- 
dent of /, albeit with a large constant prefactor, or at O (/log/) operations. A brute-force 
computation of the monodromy matrix by repeated integrations of the e.o.m would have 
complexity O (f^) because / passes are required in order to obtain / linearly independent 
solutions. Approximations have been developed which yield a fundamental matrix solution 
in O (f^)- From the matrizant the characteristic exponents can be obtained. In and 
similarly in [5,], following earlier works, a solution to recursion relations (similar to, e.g., 
eq. ( pTl) ) is achieved by truncating a linear eigenvalue problem. The convergence of the 
expansion using continued inversions is expected to be much better (this is certainly true 
for a single degree of freedom). In jsl, an expansion using Chebyshev polynomials is used 
to approximate a fundamental matrix solution. 



Appendix B: Solution of the Double-Cosine System 

As discussed in Sec. llVt the expansion in continued inversions can be extended to treat 
coupled systems of Hill equations, and here we consider the double-cosine system 



u+[A- 2Q2 cos 2t - 2Q4 cos At] u = 0. 
By substitution of the solution ansatz, eq. f lT7|) . we get the identity 



R2nC2n — Q2 (^C2n-2 + C2n+2j — Qa {C2n-A + C2n+A^ — 0, 

which gives the two recursion relations, 

QiC2n-A = —Q2C2n-2 + R2nC2n — Q2C2n+2 — Q4C'2n+4, 

and 

Q2C2n+2 = —QiC2n+A + -R2nC*2n " Q2C2n~2 — Q4C*2n-4- 



(Bl) 



(B2) 



(B3) 



(B4) 



We do not obtain the expansion to the general order, but suffice with assigning n = 1, 2, 3 
in eq. (]B3I) and n = 0,-1, —2 in eq. (1B4|) to get two expression in a form which is a 
generalization of eqs. (123|) -(124 



C2 — T2^pQ2,I3Cq, Q-2,I3C2 — Tq^/^Cq 



(B5) 



with 



-2,13 



R2 ~ Q2R4R6 ~ QaRq ^Qi ~ QaRq ^Q2RiR6 ~ Q aR~2Q A 

Q2,l3 = Q2 + Q2R4Q4 + QaRq ^Q2R4:Q4 + Q4-R-2-R-4 



(B6) 
(B7) 
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and 



Q2 + Q2R-2Q4 + Q4R-4Q2R-2Q4: + QaRaRq 



Q~2,I3 

and we have defined the matrices 

R4 = [-R4 — Q2R6 ^^2] 
R-2 = [R-2 — Q2R-IQ2] 



-1 



-Re = [Q2 + Q2R6 ^Qi] , 
R-4 = [Q2 + Q2R-IQ4] 



Tq^P — Rq — Q2R-2R-A ~ QiR-iQi — (54-^-2^2-^-2-^-4 ~ Q^R-iQi 5 (BS) 

(B9) 

(BIO) 
(Bll) 

(B12) 

which is the generahzation of eq. fl25l) . and in fact reduces to it for Q4 = 0. The mode 
vector C2 follows from eq. (IBSp . and the other three vectors in this finite expansion can be 
obtained from 

Ca = Ra [QaCq + RqC2^ , (B13) 
C.2 = R-2 [R-aCo + QiC2^ , (B14) 

C_4 = [i?-4]"' fQ2^-2 + Q4C0) . (B15) 



The characteristic exponents are obtained as zeros of the determinant of 

^2,13 = ^0,/? ~ Q-2,/3^2,/3Q2,/3, 



Appendix C: Proof of Several Results from Sec. IIIII 

In this appendix we prove some results used in Sec. ITTTl Let us first obtain an explicit 
expression for the inverse of the Floquet-Lyapunov transformation. The matrizant $ (t) of 
eq. ( 132|1 satisfies the identity 

$ (t)t j$ (t) = J, (CI) 

as can be shown by differentiating and using eq. fl32|) (the matrizant of eq. f p5|) satisfies 
a similar identity with J replaced by K). Multiplying eq. fICip on the left with J and 
rearranging a little we find 

$(t)-i = -J$ (t)t J. 

From eq. ( ITOl) we have that 

r (t)-^ = e^*r-i (0) (t) . 

Using these two expressions and substituting eq. f lTOj) we get 

r (t)-^ = -e^*r-i (0) jr-t (o) e-^*rt (t) j, (C2) 

where T~'^ = [F"^]^ as was defined already above. By substituting t — )■ t + T in the above 
equation and using the periodicity of F and the fact the B is diagonal, we find that 

[F-i(O) JF-t(0),B] =0 (C3) 
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where [, ] is the matrix commutator. We can find the exphcit form of F ^ (0) using the fact 
that U (0) is real and V (0) is purely imaginary, 

^ ' 2 V (0) (0) 
so, 

r-'(o).r-.(o). m 

with 

M = \ (0) (0)] . (C5) 

Expanding eq. f ICSP in block form, using eq. ( 126|) . we have, writing {, }^ for the matrix 
anti-commutator, 

[M + M\B]=0, {M-M\B}^ = 0, (C6) 

which together imply 

[M, 5] = 0. (C7) 

Eq. ( 1C7|) means that the invariant subspaces of M and 5 are identical, and since B is 
diagonal, we can conclude that M must be diagonal too (or, if B has degenerate eigenfre- 
quencies, M can be made diagonal by an appropriate choice of eigenvectors). We can make 
M a scalar matrix, by demanding a proper normalization of U and V. As noted at the end 
of Sec. Ill C\ the choice in eq. ( l3T|) is unique up to a matrix that commutes with B, which 
amounts to the arbitrariness of the normalization of each of the columns of U, i.e. the fact 
that each of the vectors Cq./s^ was determined only up to a constant. This allows to to rescale 
and U and V as in eq. f p9|) and obtain eq. fj^Tj) and eq. fH2|) . 

Let us find the transformed Hamiltonian of eq. fl36|) . 



n' (c, = n{u (c, , p (C, 0) + (C, 0,0- (C8) 

By use of eq. (I42|) it follows that 

r*jr = f °, n , (C9) 



and by derivating this equation we get 

r*jr + r*jr = o. (cio) 

Substituting = Fx in eq. flMl) we have 

H = ^x*r*jnrx = ^x*r*J (eb + e) x = ^x*Hx - ^x't'jtx. (cu) 

where eqs. (1C9P - (1C10P have been used to get the second line, and H is given by eq. (137|) . We 
write explicitly 

f P Q\ ^ f VV-WV V'U* - WV* \ 
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By applying dt to eq. (H5|) and substituting u = + U*(, we may write after some rear- 
ranging 

J" = Ix'^X (C13) 

where 

._(U'LU U'LU* -2iU'U-' \ 

V u^Lu Wlu* + liim-' - tWij-' J ^ ' 

and 

L = if-'V' + U-'V\ (C15) 

and in these expressions, it is understood that f/~* = dt{U~^), i.e. the time-derivative is 
apphed after matrix inversion. From eq. fH6|) we obtain the identity 

- U^U-'V' - U^L = -iU-\ (C16) 

and by derivating UU^^ = 1/ we get in addition 

f/*f7-* = -U'U-K (C17) 

Using these identities and eq. (H2l) . A is simphfied to 



Since ^f7~^t/* j = [/^f/^*, the lower-right block of the second term above is antisymmetric 

and therefore the coefficient of C*C will equal 0. By using eqs. (1C18I) . ( 1C13P and ( ICllI) in 
eq. ( ICSP we get finally 

= ^x*Hx + 3^ (t) (C19) 

with 

y = {Cw^ - ^'wX) , w = u-^dtU. (C20) 

Classically, since C and ^ commute, 3^ = 0. 

We now evaluate the integral on the r.h.s of eq. fl65|) . 



j = j d^i'f{u\c){cmo- (C21) 

Using eq. (1C11I) we can write H as a function of ^ and C, and normal-order it (i.e. put all 
^'s on the left and C's on the right) to get 

H = i'B(- : ^x*r*Jrx : -^tr {R - B} (C22) 



where R appears in the lower-left block of eq. (1C12I) . 

For any operator O we have 



(CI : O (e, C) : 10 = (CI O iC, 10 = O iC, e^'*<. 
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and we use the following result 
The first term in eq. flC22p gives 



■ n 



(C23) 



dfi'f {u I C') (CI i'BC 10 = J2 (^0^'"^) ^^'^^^ = [C'U^U-'BC - u'U-'BC] {u I C) • (C24) 

jk 



We next evaluate the upper-left block of the second term in eq. (lC22p . 



(C25) 



jk 



By the definition of the trace and the identity in eq. fllB]) we have for the second term above, 

jk jk 

= ^tr|f/t[/-* (v'U-U'Vy^ = ^tr|t/V-0f7-if/-^f7} = hr{R-iW} (C26) 
By using the formal identity resulting from eq. f H3|) . 



iT 
3^ , 



(C27) 



we can immediately see that the first term in eq. (]C25|) and the remaining terms eq. ( ]C22p 
integrate to give simply —T, and together with C = —iBC^ we can collect all the terms to 
get that 

- dtT {B -iW} 



J 



{u\0- 



(C28) 



We now show that eq. fl68|) results from the form of \N'\ as given in eq. fl63|) . This follows 
by use of trace properties and the definition of W from eq. (1C19P in the following identity, 



dt det A = ii{ (det A) A'^dtA] 



to get 



1 det (f/f/t) 



-5/4 



W\ 4det(f/f/t)-^/^ 
The solution of eq. f l69|) is given by 



dt[-\ arg det u\ = dt f -^Jmlog det U 



jj^dt det {UU^) = --tr <^ U-^U + U'^W 



--^e{tTW}. (C29) 



1^ dtdetU 



1^ tr{(detf/)f/-i9tf/} 

-2^"^ d^uT 



-Jm{tTW} (C30) 
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We now consider the integral in eq. (177j) . Using eq. (159|) we get 

(C31) 

Changing integration variables by x = U~*u and using the fact that U^^U* = U'^U^* since 
the latter is symmetric, we get 

5n,n' = ctcH' |det [/* | 6^ ("^ "'^^ O'^^^,^, , (C32) 
where X^^^ is the integral 

Tn,n' = y"rf^fexp|-if*Cx|(/7?(x))*/7$;(x) = ■ ■ ■ n^! (27r)^/2 (det 

(C33) 

and the last equality follows from eq. 12.9(1) of [26j, using the fact that C* = C^^, which 
identifies the dual polynomials of Hn (x) (denoted Gn (x) in (26|), with their complex conju- 
gates, Hn (x)*. Finally, since det C = detUy detU, and using eq. ( 163|1 . all time-dependent 
terms cancel and eq. ( ITHj) results. 
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